#include "dna.h"
#include "kmer.h"
#include "getargs.h"
#include "libheader.h"

struct PosRange {
  int start;
  int end;
};
bool comparePosRange (PosRange i,PosRange j) { return (i.start < j.start); }

struct IndexStru{
  int dna_no;
  short position;
  unsigned long code;
  IndexStru* next;
};

void delHashTable(KMer** hash_table, uint64_t hash_table_size);
int countDeg(KMer** hash_table, uint64_t hash_table_size);
int countBranchNode(KMer** hash_table, uint64_t hash_table_size);
int extendContig(KMer* tmp_kmer, vector<Contig*> &contig_vec, KMer** hash_table, uint64_t hash_table_size, int k_length);
Contig* alignRead(DNA* read, KMer** hash_table, const uint64_t hash_table_size, const int k_length);
KMer* searchKMer(const uint64_t code, KMer** hash_table, const uint64_t hash_table_size, const int k_length);
//int mergeContig(Contig* startContig, Contig* endContig, vector<Contig*> &contig_vec, KMer** hash_table, const uint64_t hash_table_size, const int k_length);
int buildHashTable(const char* dnaSeq, const int dna_length, const int k_value, KMer** hash_table, uint64_t hash_table_size);
void ConnectKMer(const int k_value, KMer** hash_table, uint64_t hash_table_size);
void extendAllContig(vector<Contig*> &contig_vec, const int k_value, KMer** hash_table, uint64_t hash_table_size);
int local_assemble(Contig* contig, vector<Contig*> &tmp_contig_vec, vector<DNA*> &read_vec, KMer** tmp_hash_table, const uint64_t hash_table_size, const int k_value);
int removeLowThres(KMer** hash_table, const uint64_t hash_table_size, int thres_m);
void removeAllBubble(const int k_value, KMer** hash_table, const uint64_t hash_table_size);
int removeBubble(uint64_t code, KMer* tmp_kmer, KMer** hash_table, const uint64_t hash_table_size, const int k_value);
double extendContigNext(uint64_t code, string &contigStr, KMer** hash_table, uint64_t hash_table_size, const int k_value, vector<KMer*> &kmer_vec);
int alignRead2(DNA* read, Contig* contig, int k_value);

int addAdjacent(KMer* kmer_ptr,KMer** hash_table, uint64_t hash_table_size, int k_length);

void print_align_pattern(const char* dna_seq, int start_pos, int print_neighbor_length, int k_length, int seq_length);
unsigned long encode1(const unsigned short* seq, int k_length);
unsigned long encode2(const unsigned short* seq, int k_length);
int getRead(vector<DNA*> &read_vec, const char* in_filepath);
//get file name in a directory
int getdir (string dir, vector<string> &files);

const char nucleotideArr[]={'A','C','G','T'};
const int nucl_num=4;

int main(int argc,const char** argv){
  time_t start_time = time(0);  
  srand(time(0));
  
  struct ParaFormat para_format[] =
  {
    { 0, "read_filepath1", "input read file", "t", 1, 's'},
    { 1, "read_filepath2", "input read file", "f", 1, 's'},
    { 2, "out_filepath", "output file name", "O", 1, 's'},
    { 3, "contig_filepath", "idba-ud contig file name", "c", 1, 's'},
    { 4, "read_length", "lenght of read (default 75)", "r", 1, 'i'},
    { 5, "insert_dist", "insert_distance (default 200)", "i", 1, 'i'},
  };
  
  Getargs args_list(argc,argv,para_format,6);
  
  if(args_list.isError()){
    args_list.print_list(cout);
    return -1;
  }
  
  const char* out_filepath = NULL;
  const char* read_filepath1 = NULL;
  const char* read_filepath2 = NULL;
  //const char* single_filepath = NULL;
  const char* contig_filepath = NULL;
  const int k_min = 21;
  const int k_max = 21;
  const int thres_m = 2;
  int insert_dist=200;
  int read_length=75;
  const int cov_range=0;
  
/*****************************************************************************/
/*                                                                           */
/* Read in parameters                                                        */
/*                                                                           */
/*****************************************************************************/
  for(int i=0;i<(int)args_list.getParaValSize();i++){ 
    switch((args_list.getPara(i))->getIndex()){
      
      //case 0: input file path
      case 0:
        read_filepath1=(args_list.getPara(i))->getStrVal();
        break;
      
      //case 1: input file path
      case 1:
        read_filepath2=(args_list.getPara(i))->getStrVal();
        break;
/*      
      //case 2: single file path
      case 2:
        single_filepath=(args_list.getPara(i))->getStrVal();
        break;
*/        
      //case 2: output file path
      case 2:
        out_filepath=(args_list.getPara(i))->getStrVal();
        break;
      
      //case 3: contig file path
      case 3:
        contig_filepath=(args_list.getPara(i))->getStrVal();
        break;
      
      //case 4: read_length
      case 4:
        read_length=(args_list.getPara(i))->getIntVal();
        break;
      
      //case 5: insert_dist
      case 5:
        insert_dist=(args_list.getPara(i))->getIntVal();
        break;
      
    }    
  }
  
  if(read_filepath1==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the read sequences (fasta) by the -t command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }
  
  if(read_filepath2==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the read sequences (fasta) by the -f command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }
  
  if(out_filepath==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the output file (txt) by the -O command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }
  
  if(contig_filepath==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the input conitg (fasta) by the -c command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }
  
/*****************************************************************************/
/*                                                                           */
/* Read in read sequences                                                    */
/*                                                                           */
/*****************************************************************************/

  vector<DNA*> in_read_vec1;
  vector<DNA*> in_read_vec2;
  getRead(in_read_vec1, read_filepath1);
  getRead(in_read_vec2, read_filepath2);
  cout<<"number of reads in 1st reads = " <<in_read_vec1.size()<<endl;
  cout<<"number of reads in 2nd reads = " <<in_read_vec2.size()<<endl;
     
  const uint64_t hash_table_size = 1073741824; //1G
  const uint64_t tmp_hash_size = 8388608; //8M         
  KMer** hash_table=new KMer*[hash_table_size];
  KMer** tmp_hash_table=new KMer*[tmp_hash_size];
  memset(hash_table, 0, hash_table_size*sizeof(KMer*));
  memset(tmp_hash_table, 0, tmp_hash_size*sizeof(KMer*));
  
  //merge two vectors
  vector<DNA*> in_read_vec;
  for(int dna_no=0;dna_no<(int)in_read_vec1.size();dna_no++){
    in_read_vec.push_back(in_read_vec1.at(dna_no));
    in_read_vec.push_back(in_read_vec2.at(dna_no));
  }
  //in_read_vec1.clear();
  //in_read_vec2.clear();

/*****************************************************************************/
/*                                                                           */
/* Assembling                                                                */
/*                                                                           */
/*****************************************************************************/
  
  vector<Contig*> contig_vec;
  
  //for each k_value
  for(int k_value=k_min;k_value<=k_max;k_value+=2){

//--------------------------------------------------------------------------------------    
// 1. build hashtable of kmer                                                          |
//-------------------------------------------------------------------------------------- 
    cout<<"k = "<<k_value<<endl;
    
    for(int dna_no=0;dna_no<(int)in_read_vec.size();dna_no++){   
      const char* dnaSeq=in_read_vec[dna_no]->getSeq();
      int dna_length=in_read_vec[dna_no]->getSeqLength();
      buildHashTable(dnaSeq, dna_length, k_value, hash_table, hash_table_size);
    }

/*    
    for(int contig_no=0;contig_no<(int)contig_vec.size();contig_no++){
      const char* dnaSeq=contig_vec[contig_no]->getSeq();
      int dna_length=contig_vec[contig_no]->getSeqLength();
      buildHashTable(dnaSeq, dna_length, k_value, hash_table, hash_table_size);
      delete contig_vec[contig_no];
    }
    
    contig_vec.clear();
*/    
    if(k_value==k_min){
      removeLowThres(hash_table, hash_table_size, thres_m);
    }
    
//--------------------------------------------------------------------------------------    
// 2. add adjacent information                                                         |
//--------------------------------------------------------------------------------------
    ConnectKMer(k_value, hash_table, hash_table_size);
    cout<<"finished building hash table"<<endl;
    
//--------------------------------------------------------------------------------------    
// 3. remove bubble                                                                    |
//--------------------------------------------------------------------------------------    
    removeAllBubble(k_value, hash_table, hash_table_size);
    cout<<"finished removing bubble"<<endl;

/*
//--------------------------------------------------------------------------------------   
// 4. find single path                                                                 |
//--------------------------------------------------------------------------------------      
    extendAllContig(contig_vec, k_value, hash_table, hash_table_size);
    cout<<"finished searching single path"<<endl;

//--------------------------------------------------------------------------------------    
// 5. local assembling                                                                 |
//-------------------------------------------------------------------------------------- 
    
    cout<<"start local assembling"<<endl;
        
    //align paired end reads
    for(int dna_no=0;dna_no<(int)in_read_vec.size();dna_no+=2){
      if(dna_no%10000==0){
        cout<<dna_no<<" out of "<<in_read_vec.size()<<endl;
      }
      //cout<<dna_no<<endl;
      Contig* startContig = alignRead(in_read_vec.at(dna_no), hash_table, hash_table_size, k_value);
      if(startContig!=NULL){
        startContig->addNeigRead(dna_no+1);
      }
      Contig* endContig = alignRead(in_read_vec.at(dna_no+1), hash_table, hash_table_size, k_value);
      if(endContig!=NULL){
        endContig->addNeigRead(dna_no);
      }//startContig->addNeigContig(endContig);
      //endContig->addNeigContig(startContig);
    }
    
    vector<Contig*> tmp_contig_vec;
    for(int i=0;i<(int)contig_vec.size();i++){
      if(i%10000==0){
        cout<<i<<" out of "<<contig_vec.size()<<endl;
        cout<<"tmp_contig_vec.size() = "<<tmp_contig_vec.size()<<endl;
      }
      if(local_assemble(contig_vec[i], tmp_contig_vec, in_read_vec, tmp_hash_table, tmp_hash_size, k_value)==0){
        contig_vec[i]=NULL;
      }
      delHashTable(tmp_hash_table, tmp_hash_size);
    }

//--------------------------------------------------------------------------------------    
// 6. Remove tigs                                                                      |
//--------------------------------------------------------------------------------------
*/    
//--------------------------------------------------------------------------------------   
// 7. Breaking contigs                                                                 |
//--------------------------------------------------------------------------------------   
    
    //reference to idna contig
    vector<DNA*> contig_dna_vec;
    getRead(contig_dna_vec, contig_filepath);
    vector<Contig*> tmp_contig_vec;
    for(int i=0;i<(int)contig_dna_vec.size();i++){
      Contig* contig=new Contig();
      contig->setSeq(contig_dna_vec[i]->getSeq());
      tmp_contig_vec.push_back(contig);
    }
    
    cout<<"finished read in idba contig"<<endl;
    
    /*
    for(int i=0;i<(int)contig_vec.size();i++){
      if(contig_vec[i]!=NULL){
        delete contig_vec[i];
      }
    }
    contig_vec.clear();
    */
    
    //for each contig, set k_mer;
    for(int i=0;i<(int)tmp_contig_vec.size();i++){
      const char* dnaSeq=tmp_contig_vec[i]->getSeq();
      int dna_length=tmp_contig_vec[i]->getSeqLength();
      
      for(int j=0;j<dna_length-k_value+1;j++){
        uint64_t code=KMer::PattoCode(&dnaSeq[j], k_value);
        KMer* tmp_kmer=searchKMer(code, hash_table, hash_table_size, k_value);
        if(tmp_kmer!=NULL){
          tmp_kmer->setContig(tmp_contig_vec[i]);
        }
      }
    }
    
    cout<<"finished setting k_mer"<<endl;
    
    //for each reads, aligned to contigs
    for(int dna_no=0;dna_no<(int)in_read_vec.size();dna_no+=2){     
      //cout<<dna_no<<endl;
      Contig* startContig = alignRead(in_read_vec.at(dna_no), hash_table, hash_table_size, k_value);
      Contig* endContig = alignRead(in_read_vec.at(dna_no+1), hash_table, hash_table_size, k_value);
      if(startContig == endContig && startContig != NULL){        
        startContig->addNeigRead(dna_no);
      }
    }
    
    cout<<"finished aligning reads"<<endl;
    
    for(int i=0;i<(int)(int)tmp_contig_vec.size();i++){
      cout<<"contig no (i) = "<<i<<endl;
      
      vector<PosRange> pos_range_vec;
      vector<int>* read_vec=tmp_contig_vec[i]->getNeigVec();
      
      //align paired end read
      for(int j=0;j<(int)read_vec->size();j+=2){
        //cout<<"i = "<<i<<" j = "<<j<<endl;
        int start=alignRead2(in_read_vec[read_vec->at(j)], tmp_contig_vec[i], k_value);
        int end=alignRead2(in_read_vec[read_vec->at(j)+1], tmp_contig_vec[i], k_value);
        //cout<<"start = "<<start<<" end = "<<end<<endl;
        bool isValid=false;
        if(start <= 0 && end >= 0){
          if(-start > end){
            int tmp_end=-start;
            start=end;
            end=tmp_end;
            isValid=true;
          }
        } else{
          if(start >= 0 && end <= 0){
            end=-end;
            if(end > start){
              isValid=true;
            }
          }
        }
        //cout<<"new start = "<<start<<" new end = "<<end<<" dist = "<<end-start+read_length<<endl;
        PosRange tmp_range;
        tmp_range.start=start;
        tmp_range.end=end;
        pos_range_vec.push_back(tmp_range);        
      }
                 
      //break contig
      if(pos_range_vec.size()<=1){
        contig_vec.push_back(tmp_contig_vec[i]);
      } else {
        sort (pos_range_vec.begin(), pos_range_vec.end(), comparePosRange);
        for(int k=0;k<(int)pos_range_vec.size();k++){
          cout<<"   ("<<pos_range_vec[k].start<<","<<pos_range_vec[k].end<<")"<<endl;
        }
        
        //determine potential break points
        vector<int> break_vec;
        const char* dnaSeq=tmp_contig_vec[i]->getSeq();
        int dna_length=tmp_contig_vec[i]->getSeqLength();
        int count_arr[dna_length-k_value+1];
        double cov_arr[dna_length-k_value+1];
      
        for(int j=0;j<dna_length-k_value+1;j++){
          uint64_t code=KMer::PattoCode(&dnaSeq[j], k_value);
          KMer* tmp_kmer=searchKMer(code, hash_table, hash_table_size, k_value);
          if(tmp_kmer!=NULL){
            count_arr[j]=tmp_kmer->getCount();
          } else {
            if(j==0){
              count_arr[j]=0;
            } else {
              count_arr[j]=count_arr[j-1];
            }
          }
        }
        
        cout<<"coverage : ";
        for(int j=0;j<dna_length-k_value+1;j++){
          int count=0;
          double sum=0;
          for(int k=max(0,j-cov_range);k<=min(j+cov_range,dna_length-k_value);k++){
            sum+=count_arr[k];
            count++;
          }
          cov_arr[j]=sum/count;
          cout<<cov_arr[j]<<",";
        }
        cout<<endl;
        
        cout<<"break points : ";
        for(int j=1;j<dna_length-k_value+1;j++){
          if(cov_arr[j] >= (double)cov_arr[j-1]*1.2 || (double)cov_arr[j]*1.2 <= cov_arr[j-1]){
            break_vec.push_back(j);
            cout<<j<<",";
          }
        }	
        cout<<endl;
        
        //deteremine repeat
        int tmp_contig_len=tmp_contig_vec[i]->getSeqLength();
        int break_num=break_vec.size();
        bool is_first_end_arr[tmp_contig_len];
        bool is_second_end_arr[tmp_contig_len];
        memset (is_first_end_arr,0,tmp_contig_len*sizeof(bool));
        memset (is_second_end_arr,0,tmp_contig_len*sizeof(bool));
        
        for(int j=0;j<(int)pos_range_vec.size();j++){
          for(int k=0;k<break_num;k++){
            for(int m=pos_range_vec[j].start;m<min(pos_range_vec[j].start+(read_length*2/3),tmp_contig_len);m++){
              is_first_end_arr[m]=true;
            }
            
            for(int m=max(0,pos_range_vec[j].end-(read_length*2/3));m<pos_range_vec[j].end;m++){
              is_second_end_arr[m]=true;
            }            
          }
        }
        
        for(int j=0;j<=tmp_contig_len;j++){
          if(is_first_end_arr[j]){
            cout<<"*";
          } else {
            cout<<"_";
          }
        }
        cout<<endl;
        for(int j=0;j<=tmp_contig_len;j++){
          if(is_second_end_arr[j]){
            cout<<"*";
          } else {
            cout<<"_";
          }
        }
        cout<<endl;
       
        int start_repeat_pos=-1;
        int end_repeat_pos=-1;
        int start_repeat_index=-1;
        int end_repeat_index=-1;
                
        if(tmp_contig_len > 1.5*insert_dist){
          int start_pos=0;
          bool isbreak=false;
          for(int j=0;j<tmp_contig_len;j++){
            if(is_first_end_arr[j]){
              start_pos=j;
              break;
            }            
          }
          for(int j=start_pos+1;j<tmp_contig_len;j++){
            if(!is_first_end_arr[j]){
              start_pos=j;
              break;
            }            
          }
          
          for(int j=0;j<break_num;j++){
            if(break_vec[j]>=start_pos){
              start_repeat_index=j;
              start_repeat_pos=break_vec[j];
              isbreak=true;
              break;
            }
            
          
/*            
              for(int k=break_vec[j];k<=max_pos;k++){
                const char* start_pattern=&(tmp_contig_vec[i]->getSeq()[k]);
                uint64_t code=KMer::PattoCode(start_pattern, k_value);
                uint64_t hash_code=KMer::reverseComplement(code, k_value);
                if(code<hash_code){
                  hash_code=code;
                }
                
                KMer* start_kmer=searchKMer(hash_code, hash_table, hash_table_size, k_value);
                
                if(start_kmer!=NULL){
                  if(!start_kmer->isSinglePrev(code)){
                    start_repeat_index=j;
                    start_repeat_pos=k;
                    isbreak=true;
                    break;
                  }
                }
              }
*/
                       
          }
          for(int j=start_repeat_pos+insert_dist/2;j<tmp_contig_len;j++){
            if(is_second_end_arr[j]){
              end_repeat_index=j;
              break;
            }
          }
          for(int j=end_repeat_index;j<tmp_contig_len;j++){
            if(!is_second_end_arr[j]){
              end_repeat_index=j;
              end_repeat_pos=j;
              break;
            }
          }
          
          if(isbreak && start_repeat_pos >=0 && end_repeat_pos >=0){
            cout<<"start_pos = "<<start_pos<<endl;
            cout<<"start_repeat_index = "<<start_repeat_index<<endl;
            cout<<"end_repeat_index = "<<end_repeat_index<<endl;
            
            cout<<"start repeat pos = "<<start_repeat_pos<<endl;
            cout<<"end repeat pos = "<<end_repeat_pos<<endl;
            cout<<"diff = "<<end_repeat_pos-start_repeat_pos<<endl;
            
            
            
            Contig* contig=new Contig();
            int contig_len=end_repeat_pos-1;
            char* contig_seq=new char[contig_len+1];
            strncpy(contig_seq, &tmp_contig_vec[i]->getSeq()[0], contig_len);
            contig_seq[contig_len]=0;
            contig->setSeq(contig_seq);
            contig_vec.push_back(contig);
            
            contig=new Contig();
            contig_len=tmp_contig_vec[i]->getSeqLength()-start_repeat_pos;
            contig_seq=new char[contig_len+1];
            strncpy(contig_seq, &tmp_contig_vec[i]->getSeq()[start_repeat_pos], contig_len);
            contig_seq[contig_len]=0;
            contig->setSeq(contig_seq);
            contig_vec.push_back(contig);
          } else {
            contig_vec.push_back(tmp_contig_vec[i]);
          }
        } else {
          contig_vec.push_back(tmp_contig_vec[i]);
        }
        
/*        
        for(int k=0;k<(int)pos_range_vec.size();k++){
          int start, end;
          if(k==0){
            start=0;
          } else{
            start = pos_range_vec[k].start;
          }
          
          if(k==(int)pos_range_vec.size()-1){
            end=tmp_contig_vec[i]->getSeqLength()-1;
          } else{
            end = pos_range_vec[k+1].start-1;
          }
          
          Contig* contig=new Contig();
          int contig_len=end-start+1;
          char* contig_seq=new char[contig_len+1];
          strncpy(contig_seq, &tmp_contig_vec[i]->getSeq()[start], contig_len);
          contig_seq[contig_len]=0;
          contig->setSeq(contig_seq);
          contig_vec.push_back(contig);
        }
        delete tmp_contig_vec[i];
*/
      }
    }
    
    cout<<"finished breaking contig"<<endl;
    
    
    
    tmp_contig_vec.clear();

    
    
    
/***************************************************************************************/    
/* 8. output contigs                                                                   */
/***************************************************************************************/ 
    
/*    
    for(int i=0;i<(int)(int)tmp_contig_vec.size();i++){
      contig_vec.push_back(tmp_contig_vec[i]);
    }
    tmp_contig_vec.clear();
*/
    
/***************************************************************************************/    
/* 9. clear HashTable                                                                  */
/***************************************************************************************/   
  
    delHashTable(hash_table, hash_table_size);
  }

/*****************************************************************************/
/*                                                                           */
/* writing contigs                                                           */
/*                                                                           */
/*****************************************************************************/  
  ofstream out_file(out_filepath);
  if (!out_file) {
    cout << "output file not found : " << out_filepath<<endl;
    return 1;
  }
  
  for(int i=0;i<(int)contig_vec.size();i++){
    out_file<<">contig"<<i<<endl;
    out_file<<contig_vec[i]->getSeq()<<endl;
  }
  out_file.close();
  
  
  //delete pointers
  for(int i=0;i<(int)in_read_vec.size();i++){
    delete in_read_vec[i];
  }
  for(int i=0;i<(int)contig_vec.size();i++){
    delete contig_vec[i];
  } 
  delete hash_table;
  delete tmp_hash_table;
  
  cout<<"Total time = "<<(time(0) - start_time)/60 <<" min"<<endl; 
  
  return 0;

}

/*
int mergeContig(Contig* startContig, Contig* endContig, vector<Contig*> &contig_vec, KMer** hash_table, const uint64_t hash_table_size, const int k_length){
  vector<Contig*> path_vec;
  if(searchContig(startContig, endContig, hash_table, hash_table_size, k_length, (int)(insert_dist*1.1), path_vec)==1){
    string contig(startContig);
    uint64_t code=startContig.getEndCode();
    for(int i=0;i<(int)path_vec.size();i++){
      if(path_vec[i].getstartCode==code){
        contig.append(&path_vec[i].getPattern()[k_length]);
        code=path_vec[i].getEndCode();
      } else {
        if(path_vec[i].getEndCodeR==code){
          contig.append(&path_vec[i].getPatternR()[k_length]);
          code=path_vec[i].getStartCodeR();
        }
      }
    }
  }
  
  
  return 0;
}
*/


int removeLowThres(KMer** hash_table, const uint64_t hash_table_size, int thres_m){
  int count=0;
  for(uint64_t i=0;i<hash_table_size;i++){
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];
      hash_table[i]=NULL;
      
      KMer* next_kmer;
      
      while(tmp_kmer!=NULL){
        next_kmer=tmp_kmer->getNextKmer();
        if(tmp_kmer->getCount()<thres_m){
          delete tmp_kmer;
          count++;
        } else {
          tmp_kmer->setNextKmer(hash_table[i]);
          hash_table[i]=tmp_kmer;
        }
        tmp_kmer=next_kmer;        
      }
    }
  }
  return count;
}

int local_assemble(Contig* contig, vector<Contig*> &tmp_contig_vec, vector<DNA*> &read_vec, KMer** hash_table, const uint64_t hash_table_size, const int k_value){

  if(contig->getNeigVec()->size()==0){
    tmp_contig_vec.push_back(contig);
    return 0;
  }

  int initial_contig_num=(int)tmp_contig_vec.size();
/***************************************************************************************/    
/* 1. build hashtable of kmer                                                          */
/***************************************************************************************/ 
  buildHashTable(contig->getSeq(), contig->getSeqLength(), k_value, hash_table, hash_table_size);

  vector<int>* neig_read_vec=contig->getNeigVec();
  for(int i=0;i<(int)neig_read_vec->size();i++){   
    const char* dnaSeq=read_vec[neig_read_vec->at(i)]->getSeq();
    int dna_length=read_vec[neig_read_vec->at(i)]->getSeqLength();
    buildHashTable(dnaSeq, dna_length, k_value, hash_table, hash_table_size);
  }
    
/***************************************************************************************/    
/* 2. add adjacent information                                                         */
/***************************************************************************************/
  ConnectKMer(k_value, hash_table, hash_table_size);
    
/***************************************************************************************/    
/* 3. remove bubble                                                                    */
/***************************************************************************************/    
  removeAllBubble(k_value, hash_table, hash_table_size);

/***************************************************************************************/    
/* 4. find single path                                                                 */
/***************************************************************************************/      
  extendAllContig(tmp_contig_vec, k_value, hash_table, hash_table_size);
  
  return (int)tmp_contig_vec.size() - initial_contig_num;
}

void extendAllContig(vector<Contig*> &contig_vec, const int k_value, KMer** hash_table, uint64_t hash_table_size){      
  for(uint64_t i=0;i<hash_table_size;i++){
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];      
      while(tmp_kmer!=NULL){
        extendContig(tmp_kmer, contig_vec, hash_table, hash_table_size, k_value);
        tmp_kmer=tmp_kmer->getNextKmer();
      }
    }
  }
}

void ConnectKMer(const int k_value, KMer** hash_table, uint64_t hash_table_size){    
  for(uint64_t hash_value=0;hash_value<hash_table_size;hash_value++){
    if(hash_table[hash_value]!=NULL){
      KMer* cur_kmer=hash_table[hash_value];
      while(cur_kmer!=NULL){
        addAdjacent(cur_kmer, hash_table, hash_table_size, k_value);
        cur_kmer=cur_kmer->getNextKmer();
      }
    }
  }
}

void removeAllBubble(const int k_value, KMer** hash_table, const uint64_t hash_table_size){
  for(uint64_t i=0;i<hash_table_size;i++){
    if(i%(hash_table_size/10)==0){
      cout<<"remove bubble i = "<<i<<" out of "<<hash_table_size<<endl;
    }
     
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];      
      while(tmp_kmer!=NULL){
        removeBubble(tmp_kmer->getPatCode(), tmp_kmer, hash_table, hash_table_size, k_value);
        removeBubble(tmp_kmer->getPatCodeR(k_value), tmp_kmer, hash_table, hash_table_size, k_value);
        tmp_kmer=tmp_kmer->getNextKmer();
      }
    }
  }
}

int buildHashTable(const char* dnaSeq, const int dna_length, const int k_value, KMer** hash_table, uint64_t hash_table_size){
  for(int i=0;i<dna_length-k_value+1;i++){
    KMer* tmp_kmer=new KMer(NULL, &dnaSeq[i], k_value);
    uint64_t code=tmp_kmer->getPatCode();
    
    KMer* insert_pos=searchKMer(code, hash_table, hash_table_size, k_value);
    if(insert_pos==NULL){
      uint64_t hash_value=code%hash_table_size;
      tmp_kmer->setNextKmer(hash_table[hash_value]);
      hash_table[hash_value]=tmp_kmer;
    } else {
      insert_pos->addCount(1);
      delete tmp_kmer;
    }
  }
  
  return 0;
}

Contig* alignRead(DNA* read, KMer** hash_table, const uint64_t hash_table_size, const int k_length){
  int read_length=read->getSeqLength();
  const char* readSeq=read->getSeq();
  vector<Contig*> contig_vec;
  vector<int> count_vec;
  
  for(int i=0;i<read_length-k_length+1;i++){
    uint64_t code=KMer::PattoCode(&readSeq[i], k_length);
    
    KMer* kmerPos=searchKMer(code, hash_table, hash_table_size, k_length);
    if(kmerPos!=NULL && kmerPos->getContig()!=NULL){
      bool isExist=false;
      for(int j=0;j<(int)count_vec.size();j++){
        if(kmerPos->getContig()==contig_vec[j]){
          isExist=true;
          count_vec[j]=count_vec[j]+1;
          break;
        }
      }
      if(!isExist){
        contig_vec.push_back(kmerPos->getContig());
        count_vec.push_back(1);
      }
    }
  }
  
  Contig* max_contig=NULL;
  int max_count=0;
  
  for(int i=0;i<(int)count_vec.size();i++){
    if(count_vec[i]>max_count){
      max_count=count_vec[i];
      max_contig=contig_vec[i];
    }
    
  }
  
  return max_contig;
}

int alignRead2(DNA* read, Contig* contig, int k_value){
  const char* read_seq=read->getSeq();
  const char* read_seqR=read->getSeqR();
  const int read_len=read->getSeqLength();
  
  const char* contig_seq=contig->getSeq();
  const int contig_len=contig->getSeqLength();
  
  int align_pos=0;
  bool is_find=false;
  
  for(int i=0;i<read_len-k_value+1;i++){
    if(is_find){
      break;
    }
    for(int j=0;j<contig_len-k_value+1;j++){
      if(strncmp(&read_seq[i], &contig_seq[j], k_value)==0){
        align_pos=j-i;
        if(align_pos<0){
          align_pos=0;
        } else {
          if(align_pos>=contig_len){
            align_pos=contig_len-1;
          }
        }
        is_find=true;
        break;
      }
    }
  }
  
  if(!is_find){
    for(int i=read_len-k_value;i>=0;i--){
      if(is_find){
        break;
      }
      for(int j=contig_len-k_value;j>=0;j--){
        if(strncmp(&read_seqR[i], &contig_seq[j], k_value)==0){
          align_pos=j-i+read_len-1;
          
          if(align_pos<0){
            align_pos=0;
          } else {
            if(align_pos>=contig_len){
              align_pos=contig_len-1;
            }
          }
          align_pos=-align_pos;
          is_find=true;
          break;
        }
      }
    }
  }
  return align_pos;
}

int removeBubble(uint64_t code, KMer* tmp_kmer, KMer** hash_table, const uint64_t hash_table_size, const int k_value){
  if(tmp_kmer->getOutDeg(code)==2){
    uint64_t code1;
    uint64_t code2;
    
    char pattern[k_value+1];
    tmp_kmer->getDNA(pattern, k_value);
    string contigStr1(pattern);    
    string contigStr2(pattern);    
    int nextNucl1=0;
    int nextNucl2=0;
        
    for(int i=0;i<4;i++){
      if(tmp_kmer->checkNext(code, i)){
        contigStr1.append(1,nucleotideArr[i]);
        nextNucl1=i;
        break;
      }
    }
    for(int i=nextNucl1+1;i<4;i++){
      if(tmp_kmer->checkNext(code, i)){
        contigStr2.append(1,nucleotideArr[i]);
        nextNucl2=i;
        break;
      }
    }
    
    code1=KMer::getNext(code, nextNucl1, k_value);
    code2=KMer::getNext(code, nextNucl2, k_value);
    vector<KMer*> kmer_vec1;
    vector<KMer*> kmer_vec2;
    
    double coverage1=extendContigNext(code1, contigStr1, hash_table, hash_table_size, k_value, kmer_vec1);
    double coverage2=extendContigNext(code2, contigStr2, hash_table, hash_table_size, k_value, kmer_vec2);
    
    int contig_len1=contigStr1.length();
    int contig_len2=contigStr2.length();
    
    if(contig_len1 == contig_len2 && strcmp(&contigStr1.c_str()[contig_len1-k_value],&contigStr2.c_str()[contig_len2-k_value])==0){
      int count_diff=0;
      for(int i=0;i<contig_len1;i++){
        if(contigStr1[i]!=contigStr2[i]){
          count_diff++;
        }
      }
      if(count_diff<2){
      //delete one of the path
        cout<<"mergining"<<endl;
        cout<<contigStr1<<endl;
        cout<<contigStr2<<endl;
        
        if(coverage1 < coverage2){
          tmp_kmer->setNext(code, nextNucl1, false);
          uint64_t last_code=KMer::PattoCode(&contigStr1.c_str()[contig_len1-k_value], k_value);
          kmer_vec1[kmer_vec1.size()-1]->setPrev(last_code,contigStr1.c_str()[contig_len1-k_value-1],false);
          
          for(int i=0;i<(int)kmer_vec1.size()-1;i++){
            uint64_t hash_code=kmer_vec1[i]->getPatCode();
            KMer* prev_kmer=hash_table[hash_code%hash_table_size];
            if(hash_code==hash_table[hash_code%hash_table_size]->getPatCode()){
              hash_table[hash_code%hash_table_size]=prev_kmer->getNextKmer();
              delete prev_kmer;
            } else {
              KMer* cur_kmer=prev_kmer->getNextKmer();
              while(cur_kmer!=NULL && cur_kmer->getPatCode()!=hash_code){
                prev_kmer=cur_kmer;
                cur_kmer=cur_kmer->getNextKmer();
              }
              if(cur_kmer!=NULL){
                prev_kmer->setNextKmer(cur_kmer->getNextKmer());
                delete cur_kmer;
              }
            }
          }          
        } else {
          tmp_kmer->setNext(code, nextNucl2, false);
          uint64_t last_code=KMer::PattoCode(&contigStr2.c_str()[contigStr2.length()-k_value], k_value);
          kmer_vec2[kmer_vec2.size()-1]->setPrev(last_code,contigStr2.c_str()[contig_len2-k_value-1],false);
          
          for(int i=0;i<(int)kmer_vec2.size()-1;i++){
            uint64_t hash_code=kmer_vec2[i]->getPatCode();
            KMer* prev_kmer=hash_table[hash_code%hash_table_size];
            if(hash_code==hash_table[hash_code%hash_table_size]->getPatCode()){
              hash_table[hash_code%hash_table_size]=prev_kmer->getNextKmer();
              delete prev_kmer;
            } else {
              KMer* cur_kmer=prev_kmer->getNextKmer();
              while(cur_kmer!=NULL && cur_kmer->getPatCode()!=hash_code){
                prev_kmer=cur_kmer;
                cur_kmer=cur_kmer->getNextKmer();
              }
              if(cur_kmer!=NULL){
                prev_kmer->setNextKmer(cur_kmer->getNextKmer());
                delete cur_kmer;
              }
            }
          }          
        }
      }
    }
    return 1;
  } else {
    return 0;
  }
}

double extendContigNext(uint64_t code, string &contigStr, KMer** hash_table, uint64_t hash_table_size, const int k_value, vector<KMer*> &kmer_vec){
  uint64_t hash_code=KMer::reverseComplement(code, k_value);
  if(code<hash_code){
    hash_code=code;
  }
  
  KMer* cur_kmer=searchKMer(hash_code, hash_table, hash_table_size, k_value);
  if(cur_kmer!=NULL){
    int countNum=cur_kmer->getCount();
    
    kmer_vec.push_back(cur_kmer);
    Contig* dummyContig=new Contig();
    cur_kmer->setContig(dummyContig);
    
    bool isSingleNext=cur_kmer->isSingleNext(code);
    while(isSingleNext){
      int nextNucl=0;
      for(int i=0;i<4;i++){
        if(cur_kmer->checkNext(code, i)){
          contigStr.append(1,nucleotideArr[i]);
          nextNucl=i;
          break;
        }
      }
      
      code=KMer::getNext(code, nextNucl, k_value);
      uint64_t codeR=KMer::reverseComplement(code, k_value);
      hash_code=(code<codeR)?code:codeR;
      
      cur_kmer=searchKMer(hash_code, hash_table, hash_table_size, k_value);
      	
      if(cur_kmer->getContig()!=NULL){
        break;
      } else {
        isSingleNext=cur_kmer->isSingleNext(code);
        cur_kmer->setContig(dummyContig);
        kmer_vec.push_back(cur_kmer);
        countNum+=cur_kmer->getCount();
        if(!cur_kmer->isSinglePrev(code)){
          break;
        }
      }      
    }
    
    for(int i=0;i<(int)kmer_vec.size();i++){
      kmer_vec[i]->setContig(NULL);
    }
    delete dummyContig;
    
    return ((double)countNum)/kmer_vec.size();
  } else {
    cout<<"Null cur_kmer in extendContigNext"<<endl;
    return 0;
  }
}

KMer* searchKMer(const uint64_t code, KMer** hash_table, const uint64_t hash_table_size, const int k_length){
  KMer* kmer=hash_table[code%hash_table_size];
  while(kmer!=NULL){
    if(kmer->getPatCode()==code){
      return kmer;
    } else {
      kmer=kmer->getNextKmer();
    }
  }
  return NULL;
}

int extendContig(KMer* tmp_kmer, vector<Contig*> &contig_vec, KMer** hash_table, uint64_t hash_table_size, int k_length){
  if(tmp_kmer->getContig()!=NULL){
    return 0;
  }
  
  uint64_t code=tmp_kmer->getPatCode();
  bool isSingleNext=tmp_kmer->isSingleNext(code);
  bool isSinglePrev=tmp_kmer->isSinglePrev(code);
    
  if(isSingleNext && isSinglePrev){
    Contig* contigPtr=new Contig();
     
    char pattern[k_length+1];
    tmp_kmer->getDNA(pattern, k_length);
    string contigStr(pattern);    
    KMer* startKmer=tmp_kmer;
    tmp_kmer->setContig(contigPtr);
        
    while(isSingleNext){
      int nextNucl=0;
      for(int i=0;i<4;i++){
        if(tmp_kmer->checkNext(code, i)){
          contigStr.append(1,nucleotideArr[i]);
          nextNucl=i;
          break;
        }
      }
      
      code=KMer::getNext(code, nextNucl, k_length);
      uint64_t codeR=KMer::reverseComplement(code, k_length);
      uint64_t hash_code=(code<codeR)?code:codeR;
      
      tmp_kmer=searchKMer(hash_code, hash_table, hash_table_size, k_length);
      	
      if(tmp_kmer->getContig()!=NULL){
        break;
      } else {
        isSingleNext=tmp_kmer->isSingleNext(code);
        tmp_kmer->setContig(contigPtr);
        if(!tmp_kmer->isSinglePrev(code)){
          break;
        }
      }      
    }
    
    string preContigStr("");
    tmp_kmer=startKmer;
    code=tmp_kmer->getPatCode();
    
    while(isSinglePrev){
      int prevNucl=0;
      for(int i=0;i<4;i++){
        if(tmp_kmer->checkPrev(code, i)){
          preContigStr.append(1,nucleotideArr[i]);
          prevNucl=i;
          break;
        }
      }
      
      code=KMer::getPrev(code, prevNucl, k_length);
      uint64_t codeR=KMer::reverseComplement(code, k_length);
      uint64_t hash_code=(code<codeR)?code:codeR;
      
      tmp_kmer=searchKMer(hash_code, hash_table, hash_table_size, k_length);
      
      if(tmp_kmer->getContig()!=NULL){
        break;
      } else {
        isSinglePrev=tmp_kmer->isSinglePrev(code);
        tmp_kmer->setContig(contigPtr);
        if(!tmp_kmer->isSingleNext(code)){
          break;
        }
      }
    }
    
    int pre_length=preContigStr.length();
    int after_length=contigStr.length();
    int contig_len=pre_length + after_length;
    char* contig_pat=new char[contig_len+1];
    for(int i=0;i<pre_length;i++){
      contig_pat[i]=preContigStr.at(pre_length-1-i);
    } 
    for(int i=0;i<after_length;i++){
      contig_pat[i+pre_length]=contigStr.at(i);
    }
    contig_pat[contig_len]=0;
    
    contigPtr->setSeq(contig_pat);
    contig_vec.push_back(contigPtr);
    
    char* startPat=new char[k_length+1];
    startKmer->getDNA(startPat, k_length);
    //cout<<startPat<<" "<<contig_pat<<endl;
    delete startPat;
    
    return contig_len;
  }
  
  return 0;
}

int addAdjacent(KMer* kmer_ptr,KMer** hash_table, uint64_t hash_table_size, int k_length){
  uint64_t pat_code=kmer_ptr->getPatCode();
  for(int i=0;i<4;i++){
    uint64_t next_code = KMer::getNext(pat_code, i, k_length);
    uint64_t next_codeR = KMer::reverseComplement(next_code, k_length);
    uint64_t hash_code = next_code < next_codeR ? next_code : next_codeR;
    
    uint64_t hash_value = hash_code % hash_table_size;
        
    KMer* cur_kmer=hash_table[hash_value];
    while(cur_kmer!=NULL){
      if(cur_kmer->getPatCode()==hash_code){
        kmer_ptr->setNext(pat_code, i, true);
        break;
      } else {
        cur_kmer=cur_kmer->getNextKmer();
      }
    }
  }
  
  for(int i=0;i<4;i++){
    uint64_t prev_code = KMer::getPrev(pat_code, i, k_length);
    uint64_t prev_codeR = KMer::reverseComplement(prev_code, k_length);
    uint64_t hash_code = prev_code < prev_codeR ? prev_code : prev_codeR;
    
    uint64_t hash_value = hash_code % hash_table_size;
        
    KMer* cur_kmer=hash_table[hash_value];
    while(cur_kmer!=NULL){
      if(cur_kmer->getPatCode()==hash_code){
        kmer_ptr->setPrev(pat_code, i, true);
        break;
      } else {
        cur_kmer=cur_kmer->getNextKmer();
      }
    }
    
  }
  return 0;
}

int countBranchNode(KMer** hash_table, uint64_t hash_table_size){
  int count_branch=0;
  
  for(uint64_t i=0;i<hash_table_size;i++){
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];
     
      while(tmp_kmer!=NULL){
        int first_four_bit = (int)tmp_kmer->getAdjArr()[0]+(int)tmp_kmer->getAdjArr()[1]
                             +(int)tmp_kmer->getAdjArr()[2]+(int)tmp_kmer->getAdjArr()[3];
        if( first_four_bit >= 2){
          count_branch++;   
        }
        
        if((int)(tmp_kmer->getAdjArr().count()) - first_four_bit >= 2){
          count_branch++;
        }
        tmp_kmer=tmp_kmer->getNextKmer();
      }
    }
  }
  
  return count_branch;
}

int countDeg(KMer** hash_table, uint64_t hash_table_size){
  int count_deg=0;
  
  for(uint64_t i=0;i<hash_table_size;i++){
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];
     
      while(tmp_kmer!=NULL){
        count_deg+=(int)(tmp_kmer->getAdjArr().count());
        tmp_kmer=tmp_kmer->getNextKmer();
      }
    }
  }
  
  return count_deg;
}

void delHashTable(KMer** hash_table, uint64_t hash_table_size){
  for(uint64_t i=0;i<hash_table_size;i++){
    if(hash_table[i]!=NULL){
      KMer* tmp_kmer=hash_table[i];
      KMer* next_kmer;
      
      while(tmp_kmer!=NULL){
        next_kmer=tmp_kmer->getNextKmer();
        delete tmp_kmer;
        tmp_kmer=next_kmer;
      }
      hash_table[i]=NULL;
    }
  }
}

void print_align_pattern(const char* dna_seq, int start_pos, int print_neighbor_length, int k_length, int seq_length){
  if(start_pos>=print_neighbor_length){
    cout<<"...";
    for(int k=print_neighbor_length;k>0;k--){
      cout<<dna_seq[start_pos-k];
    }    
  } else{
    cout<<"   ";
    for(int k=0;k<print_neighbor_length-start_pos;k++){
      cout<<" ";
    }
    for(int k=0;k<start_pos;k++){
      cout<<dna_seq[start_pos-k];
    }
  }
  cout<<"   ";
  for(int k=0;k<k_length;k++){
    cout<<dna_seq[start_pos+k];
  }
  cout<<"   ";
  if(start_pos+k_length+print_neighbor_length<seq_length){
    for(int k=start_pos+k_length;k<start_pos+k_length+print_neighbor_length;k++){
      cout<<dna_seq[k];
    }
    cout<<"..."<<endl;
  } else {
    for(int k=start_pos+k_length;k<seq_length;k++){
      cout<<dna_seq[k];
    }
    cout<<endl;
  }
}

unsigned long encode1(const unsigned short* seq, int k_length){
  unsigned long value1=0x0;
  
  for(int i=0;i<k_length;i++){
    value1 = value1 << 2;
    value1 = value1 | (unsigned long)seq[i];
  }
  
  return value1;
}

unsigned long encode2(const unsigned short* seq, int k_length){
  unsigned long value2=0x0;
  
  for(int i=k_length-1;i>=0;i--){
    value2 = value2 << 2;
    value2 = value2 | (unsigned long)(0x3 - seq[i]);
  }
  
  return value2;
}

int getRead(vector<DNA*> &read_vec, const char* in_filepath){
  ifstream read_file(in_filepath);
  if (!read_file) {
    cout << "prot file not found : " << in_filepath<<endl;
    return 1;
  }
   
  const int buffer_len = 67108864;
  char* line;
  char* line_buf=new char[buffer_len];
  char* seq_buf=new char[buffer_len];
  char* seq_name=new char[buffer_len];
  seq_name[0]=0;
  string tmp_seq("");
  
  DNA *dna = NULL;
  while (!read_file.eof()) {
    read_file.getline(line_buf,buffer_len);
    
    int start=0;
    while(isspace(line_buf[start]) && start < (int)strlen(line_buf)){
      start++;
    }
    
    if(start==(int)strlen(line_buf))
      continue;
    else
      line=&line_buf[start];
    
    if (line[0] == '>') {
      //read a new sequence
      int seq_len = tmp_seq.length();
      if ( seq_len > 0 ) {
	dna = new DNA( strdup(seq_name), strdup(tmp_seq.c_str()) );
	read_vec.push_back(dna);
	tmp_seq = "";
      }
      strncpy(seq_name, &line[1], strlen(line)-1);
    } else {
      int end = buffer_len - start;
      int current = 0;
      
      while (current < end && line[current] != 0) {
	if ( line[current] == '#' )
	  break; // ignore comment
	char c = toupper(line[current]);
	if ( 'A' <= c && c <= 'Z' ) {
	  seq_buf[current]=c;
	  current ++;
	} else {
	  cout<<"In main.cpp: incorrect character \""<<line[current]<<"\""<<endl;
	  break;
	}
	
      }
      seq_buf[current]=0;
      tmp_seq += seq_buf;
    }
  }

  // last sequence
  dna = new DNA( strdup(seq_name), strdup(tmp_seq.c_str()) );
  read_vec.push_back(dna);
  tmp_seq = "";
  //quence seq(seq_name, tmp_seq.c_str() );
  //quences.push_back(seq);
  //q_file.close();

  read_file.close();
  
  delete[] line_buf;
  delete[] seq_buf;
  delete[] seq_name;
  
  return 0;
}

//get file name in a directory
int getdir (string dir, vector<string> &files)
{
    DIR *dp;
    struct dirent *dirp;
    if((dp  = opendir(dir.c_str())) == NULL) {
        cout << "Error(" << errno << ") opening " << dir << endl;
        return errno;
    }

    while ((dirp = readdir(dp)) != NULL) {
        files.push_back(string(dirp->d_name));
    }
    closedir(dp);
    return 0;
}
